#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

#define SQUARE(x) ((x) * (x))
#define SQRT(x) (sqrt(x))
#define STUDENT 1.664           // for 80 experiments (n=79)


typedef struct stat
{
   double mean;
   double variance;
} stat_t;

typedef struct data
{
   stat_t* x;
   stat_t* y;
   stat_t* z;
   int n;
} data_t;

typedef struct point
{
   double x;
   double y;
   double z;
} point_t;

void print_usage()
{
   printf("stat fileprefix numfiles numpoints resultfile\n");
}

inline void calc_mean(int index, point_t* p, data_t* d)
{
   data_t* dp = &d[index];

   dp->x->mean += p->x / dp->n;
   dp->y->mean += p->y / dp->n;
}

inline void calc_var(int index, point_t* p, data_t* d)
{
   data_t* dp = &d[index];

   dp->x->variance += SQUARE(p->x - dp->x->mean) / dp->n;
   dp->y->variance += SQUARE(p->y - dp->y->mean) / dp->n;
}

void read_file(char* prefix, int n, int size, point_t* pp, data_t* d, int mean)
{
   FILE* f;
   char s[100];
   char num[10];
   point_t *p;
   int i;

   strncpy(s, prefix, 100);
   sprintf(num, "%d", n);
   strcat(s, num);

   f = fopen(s, "r");

   if (!f)
   {
      puts("Invalid file");
      exit(-1);
   }

   for (i = 0; i < size; i++)
   {
      p = &pp[i];

      if (fscanf(f, "%lf\t%lf\t%lf\n", &p->x, &p->y, &p->z) != 3)
      {
         puts("Invalid data");
         exit(1);
      }

      if (mean)
         calc_mean(i, p, d);
      else
         calc_var(i, p, d);
   }

   fclose(f);
}

inline double get_confidence_interval(double mean, double variance, int n)
{
   return ( ( STUDENT * SQRT(variance) ) / SQRT(n) );
}

int main(int argn, char** argv)
{
   char* fileprefix;
   int numfiles;
   int n;
   data_t* d;
   data_t* dd;
   point_t* p;
   int i;

   if (argn != 4) 
   {
      print_usage();
      return -1;
   }

   fileprefix = argv[1];
   numfiles = atoi(argv[2]);
   n = atoi(argv[3]);

   if (n <= 0)
   {
      puts("Must have more than 0 points.");
      return -1;
   }

   d = (data_t*)malloc(sizeof(data_t)*n);
   
   for (i=0; i< n; i++)
   {
      d[i].n = numfiles;
      d[i].x = (stat_t*)malloc(sizeof(stat_t)*numfiles);
      d[i].y = (stat_t*)malloc(sizeof(stat_t)*numfiles);
      bzero(d[i].x, sizeof(stat_t)*numfiles);
      bzero(d[i].y, sizeof(stat_t)*numfiles);
   }

   p = (point_t*)malloc(sizeof(point_t)*n);

   for (i=1; i<= numfiles; i++)
   {
      //mean
      read_file(fileprefix, i, n, p, d, 1);
   }

   for (i=1; i<= numfiles; i++)
   {
      // variance
      read_file(fileprefix, i, n, p, d, 0);
   }

   printf("#x_mean\ty_mean\tx_error\ty_error\n");

   for (i = 0; i < n; i++)
   {
      dd = &d[i];
      printf("%lf\t%lf\t%lf\t%lf\n", dd->x->mean,
             dd->y->mean, get_confidence_interval(dd->x->mean, dd->x->variance, dd->n), 
             get_confidence_interval(dd->y->mean, dd->y->variance, dd->n));
   }

   for (i=0; i< n; i++)
   {      
      free(d[i].x);
      free(d[i].y);
   }

   free(d);
   free(p);

   return 0;
}
